A hybrid multi objective cellular spotted hyena optimizer for wellbore trajectory optimization

Cost and safety are critical factors in the oil and gas industry for optimizing wellbore trajectory, which is a constrained and nonlinear optimization problem. In this work, the wellbore trajectory is optimized using the true measured depth, well profile energy, and torque. Numerous metaheuristic algorithms were employed to optimize these objectives by tuning 17 constrained variables, with notable drawbacks including decreased exploitation/exploration capability, local optima trapping, non-uniform distribution of non-dominated solutions, and inability to track isolated minima. The purpose of this work is to propose a modified multi-objective cellular spotted hyena algorithm (MOCSHOPSO) for optimizing true measured depth, well profile energy, and torque. To overcome the aforementioned difficulties, the modification incorporates cellular automata (CA) and particle swarm optimization (PSO). By adding CA, the SHO’s exploration phase is enhanced, and the SHO’s hunting mechanisms are modified with PSO’s velocity update property. Several geophysical and operational constraints have been utilized during trajectory optimization and data has been collected from the Gulf of Suez oil field. The proposed algorithm was compared with the standard methods (MOCPSO, MOSHO, MOCGWO) and observed significant improvements in terms of better distribution of non-dominated solutions, better-searching capability, a minimum number of isolated minima, and better Pareto optimal front. These significant improvements were validated by analysing the algorithms in terms of some statistical analysis, such as IGD, MS, SP, and ER. The proposed algorithm has obtained the lowest values in IGD, SP and ER, on the other side highest values in MS. Finally, an adaptive neighbourhood mechanism has been proposed which showed better performance than the fixed neighbourhood topology such as L5, L9, C9, C13, C21, and C25. Hopefully, this newly proposed modified algorithm will pave the way for better wellbore trajectory optimization.


Introduction
The increasing demand for energy consumption around the world has promoted the depletion of conventional energy sources. Hence, the high energy demand has drawn attention to producing energy from unconventional sources. However, producing oil/gas from unconventional sources is not easier with conventional production methods. Now, with the advancement of technologies, such as directional/horizontal drilling, the production from unconventional sources has become possible with some unwanted challenges [1]. One of the major challenges in directional / horizontal drilling is wellbore trajectory design, which is associated with cost and safety [2]. Sometimes, directional/horizontal drilling needs high expenditure, which tends to increase the oil and gas price at the consumer level. Therefore, to minimize the oil/gas price, it is crucial to minimize the operational expenditure. In directional/horizontal drilling, one of the key ways to reducing operational expenditure is optimizing the wellbore trajectory [3].
Wellbore trajectory is the direction in which the wellbore is drilled. To reach a sub-surface target, there can be thousands of probable well paths. However, the success of directional drilling depends on choosing the best path, which can be only done by the optimization of wellbore trajectory. A wellbore trajectory can be optimized considering several parameters; among them, length, torque, energy, rate of penetration, separation factor is the most influential [4][5][6][7][8][9]. Some pieces of research optimized the wellbore trajectory by considering one parameter [7]. However, this single objective optimization could not be able to provide enough cost efficiency and safety to the wellbore. Therefore, to increase cost efficiency and safety, multi-objective optimization was introduced by several researchers [5,[10][11][12][13]. In multi-objective optimization, two or more parameters were considered for wellbore trajectory optimization. Each of these parameters is optimized considering several tuning variables such as azimuth angle, dogleg severity, inclination angle, and kick-off point.
To optimize the above-mentioned parameters several algorithms were utilized with some drawbacks such as less exploitation capability, local optima trapping, nonuniformed distribution of non-dominated solutions, and disability of tracking isolated minima [11]. The researchers focused on metaheuristic algorithms due to the weakness of the traditional algorithms in the large search region [14,15]. Among the metaheuristic algorithms, genetic algorithms (GA), particle swarm optimization (PSO), ant colony optimization (ACO), artificial bee colony optimization (ABC), harmony search (HS) were utilized for well trajectory optimization [2,4,16]. To improve the issues faced by the metaheuristic algorithms and to improve the efficiency, some hybrid algorithms, for example, hybrid cuckoo search optimization (HCSO), hybrid bat flight optimization (HBFO) were introduced [17][18][19]. Due to the hybridization, these algorithms showed some significant improvements in the exploration capabilities [11]. However, these improvements enabled the algorithms to provide better solutions but on the other side, those make the convergence speed slower. Therefore, still, improvement is indispensable for increasing the exploitation capabilities of these algorithms.
In recent times, multi-objective genetic algorithm (MOGA), multi-objective cellular particle swarm optimization (MOCPSO), multi-objective cellular grey wolf optimization and particle swarm optimization (MOCGWOPSO) have been used for length, torque, and well-profile energy optimization [5,10,11]. However, MOGA and MOCPSO have faced exploitation related problems. On the other side, MOCGWOPSO provided excellent non-dominated solutions for length and torque, but it showed weakness in the case of well-profile energy optimization during multi-objective optimization. From the above discussion, it can be concluded that hybridization of any standard algorithms makes the algorithms more effective and efficient for wellbore trajectory optimization. Hence, this research focuses on the hybridization of the spotted hyena optimization (SHO) algorithm for wellbore trajectory optimization.
SHO is a newly designed metaheuristic algorithm that is inspired by the hunting mechanism of the spotted hyenas [20,21]. This algorithm has a very high convergence rate. Therefore, this algorithm faced local optima trapping issue during non-linear optimization. To overcome this issue herein cellular automata (CA) has been incorporated with SHO in this work [22]. CA is incorporated due to its slow diffusion mechanism and information exchanging capability among the neighbours. The slow diffusion mechanism helps to avoid the local optima trapping issue and the information exchanging capabilities enhance the local search capability of SHO. Besides, the velocity update mechanism of PSO has been incorporated to enhance the hunting capability of SHO. Moreover, an adaptive neighbourhood mechanism has been proposed and compared with the fixed neighbourhood topology structure such as L5, L9, C9, C13, C21, and C25 in this work. The performance analysis of the proposed algorithms for wellbore trajectory optimization has been done by comparing with the previously used MOCPSO and other states of the art algorithms such as MOSHO, MOCGWO [21,23]. This comparative analysis has been conducted based on some statistical analysis. The proposed algorithm has achieved the lowest value of IGD, SP, ER. This indicates that the obtained Pareto front by the proposed algorithm is nearer to the true Pareto front, the non-dominated solutions are nearly spaced, and it has a very less number of isolated minima. The Spearman correlation coefficient test has also been performed to analyze the sensitivity of each decision variable on the three objectives [24]. This proposed algorithm will pave the way to design a less complex and cost-effective wellbore trajectory.

Mathematical formulation
Up to the date, several methods (radius of curvature method (RCM), tangential method, angle averaging method, minimum curvature method) were utilized for characterizing the directional well design parameters [25,26]. In this work, the RCM method is used to formulate the three-objective functions [25,27].
To compute the trajectory length, radius curvature method was utilized considering several parameters such as azimuth angle, hold angle, vertical inclination, lateral length, true vertical depth, and dogleg severity. The following formulas are used to compute the constant of curvature and radius of curvature between two points in RCM.
ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi DM ¼ r � ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi Herein, a and r represents the constant of curvature and radius of curvature, respectively. The dogleg severity, inclination angle, and azimuth angle are denoted by T, ; i , and θ i , respectively. ΔM is the 3D well path between two points. The general vertical plane for a wellbore trajectory with various straight and curved sections is depicted in Fig 1. According to Fig 1, the TMD's fundamental equation is composed of seven segments which are represented by Eq (4).
Where D kop represents the estimated length of kick-of section, build and drop segments are represented by (D 1 , D 3 and D 5 ), tangent segment, hold segment and horizontal section are represented by (D 2 ), (D 4 ) and HD respectively. The torque and well-profile energy can be characterized by the following equations.
In this work buoyancy factor B = 0.7, weight of unit pipe length w = 0.3KN/ft,friction factor u = 0.2 and pipe diameter D = 0.2ft have been used for the calculation. The detail description of mathematical formulation was demonstrated in our previous work [11].

Constraints
The optimization of the wellbore trajectory is constrained in two ways: by operational constraints and by the upper and lower limits of 17 tuning variables. The tuning variables' upper and lower limits are listed in Table 1. However, in this work the operational constraints are the true vertical depth (TVD), and the casing setting depth, C1, C2, C3. During the trajectory optimization, the following equation should also be satisfied.
Herein, the vertical depth of each subsection at each drop off point is denoted by the symbol Y i2 (1,6) . There are also some non-negative constraints [10]. Fig 2 illustrates the trajectory's deviated direction to the east, north, and vertical sides. The offset distance in these three directions is derived following the radius of curvature method.

Variables
Variable constraints imposed on Wellbore design

PLOS ONE
on cost. Therefore, it is necessary to consider the following constraints during trajectory optimization.

Framework for wellbore trajectory optimization
In this multi-objective optimization, the proposed hybrid algorithm needs to optimize length, torque, and well profile energy from Eqs (4-6). The algorithm needs to fine-tune all the 17 tuning variables for optimizing these three objectives parallelly. Table A1 in S1 Appendix tabulated the bounds and explanation of these variables. The overall multi-objective optimization process can be mathematically expressed as follows.
arg min i2ð1;2;3Þ f i ðXÞ; X 2 fx j g where j ¼ 1; 2; . . . 17 In the proposed algorithm diffusion mechanism of CA will be used to increase the exploration capability of SHO, later velocity update equation of PSO will be used. This will improve the hunting mechanism of SHO.

Spotted hyena optimizer
The mechanism of SHO will be discussed in this section. It is a recently proposed metaheuristic algorithm that was inspired by the hyena's prey hunting mechanism [20,21]. Usually, the SHO algorithm imitates the behaviours of a consistent hyena cluster. It consists of four main steps such as searching, encircling, hunting and attacking. The group of fellow peers is directing the hunting behaviour toward the best solution, and preserve the best results. The following equations are used to replicate the encircling interactions of spotted hyenas.
Where M h �! define the distance to the prey from the spotted hyena, x represents the present iteration, A ! and F ! are the coefficient vectors, Q p �! specifies the prey's position vector, Q ! indicates the hyena's position vector. Besides || and "." respectively represent the absolute value and vector multiplication. The following formulas are used to calculate the vector A ! and F ! . Where On the other hand a hyena can adjust its position around its prey by applying Eqs (14)(15)(16). This algorithm stores the best solutions and compels others to upgrade their positions. When the value becomes jFj ! < 1 then the algorithm enables the hyenas to attack the target. The following equations are defined to replicate the hunting behaviour of hyenas and to identify the feasible search space regions.
Where Q h �! denotes the first best-spotted hyena's position, Q k �! defines the specific location of other search agents. The number of spotted hyenas is represented by N which is calculated as follows: Where G ! is a random vector in [0.5,1], nos denotes the number of solutions after G ! is added However in a certain search space these solutions are almost the same as the best optimum solutions. Moreover D h �! is a group of N number of optimal solutions.
Exploitation. Like every algorithm SHO also has to perform exploration and exploitation.
During exploitation, the value of h ! gradually decreases from 5 to 0. The variations in F ! also contribute to exploitation. The algorithm allows the hyenas to attack the prey when it becomes jFj ! < 1. The formulation can be expressed as follows: Where Q ! ðx þ 1Þ registers the best position and helps the search agent to update their position according to the best search agent. It permits all of its search agents to attack the prey. when the termination condition is matched, the positions of search agents which create a cluster are considered as the optimal solutions. Archive. In multi-objective optimization, the stock of all Pareto optimal solutions (those that have been obtained so far) is defined as an archive. It is comprised of two major components: an archive controller and a grid system.
Archive controller. The primary responsibility of this controller is to determine whether or not to include the solution in the archive. The following are some of the most important points to note about the updating mechanism: • If any single member of the archive dominates it, the solution cannot be included.
• If a new solution is superior to one or more solutions in the archive, it can be incorporated.
• If the new solution and archive solutions do not have any dominance over each other, the new member should be incorporated in the archive.
• The grid method should be used so that one of the crowded solution sections can remove. This will allow to include a new solution that will enhance Pareto's optimal diversity.
Grid mechanism. Grid is a space in which each particle obtains a unique solution to its objective function. An adaptive grid can be considered of as a hypercubic space with a uniform distribution of elements. If the inserted individual is outside the grid's current bounds, the grid must be reconfigured. The grid mechanism divides the search space into multiple subsearch zones.
Group selection mechanism. The challenge of comparing the solutions with archive members in a multi-objective search area is very complex. This process has been performed in this study through the use of a group selection mechanism. The mechanism determines the least populated area inside the solution region. Later in the chapter, it suggests one of the non-dominated solutions from the least populated region to the closest neighbour region. The roulette-wheel method was employed to assist in determining the least-populated area [28]. This method can be defined as where e stands for a constant number (e>1), the total number of Pareto optimal solutions achieved from the n th segment till now is represented by S, n is the number of segments This is a proportionate selection strategy. According to this method, every individual's fitness value is almost correlated to the range offered by the roulette wheel proportion. G ! is the deciding vector in MOSHO, and it is directly proportional to the number of optimal solution selections. The archive protects the obtained non-dominated solutions against degradation. The conventional algorithm, which is based on the archive, employs a variety of operators (mutation, crossover), which compelled the algorithm to focus more on the archive member. As a result, the variables in MOSHO exchange information with the various solutions in the search space. This enhances the capability for exploration but reduces the capability for convergence. That is why a group selection mechanism has been used to select a minimum of one member from the solution space. Cellular automata. The concept of cellular automata has been described in this section. Von Neumann and Ulam first published the concept of CA [11,22]. It can be characterized as a distribution of cells excreted in a particular topological structure. They are the cell, the cell state, the cell space, the neighbourhood, and the transition rule [29]. The next status of each cell will be determined by considering the present state of all the neighbouring cells. CA is composed of five components. They are cell, cell state, cell space, neighbourhood, and transition rule [30]. The cell state is a term that refers to the information about the present cell. It assists in determining the next state. Cell space is a collection of cell sets. It has applications in multiple dimensions (one, two, and three). However, because real-world processes are mostly reproduced using finite grids, the cell space boundary must be defined during operation. The limit is a ring grid. That is, the left border will remain connected to the right, and the top border will remain connected to the bottom border. The neighbourhood can be described as a group of cells that surround a centre cell. It is primarily responsible for selecting the following state. The transition rule determines the cell's next state based on the status of neighbouring cells. The formulation of CA can be characterized as follows.
Usually, an m-dimensional CA can be characterized as a grid of m-dimensional single cells. Each cell has its own value. According to the transition rule, each cell can update its state. Therefore, Q is a cellular automaton that can be formulated as Eq (23).
Where, The finite set of state :T Dimension of Q :d 2 Z + Neighborhood :H Transition rule :g Let's take i 2 Z m is the position of a cell where m is the dimension of the latticegrid, then the neighbourhood H can be formulated as Neighbourhood. The term neighbourhood refers to a group of cells that surround a central cell. Additionally, the neighbours can be defined as the other atoms connected by a single atom. The following definitions are provided to illustrate the precise composition of various neighbours. This work defines the grid's neighbours based on their direction and radius. It is frequently referred to by two structural labels. They are L n and C n . If there are n-1 nearest neighbour around the centre cell which are specific in direction, then the structure can be denoted by L n or C n .
If the directions of neighbouring cells are at the top, bottom, left, and right then it will be denoted by L n . If on the other hand, the directions are at the top, bottom, left, right, and diagonal then it will be denoted by C n . Six distinct forms of classical structure were analysed in this study to determine the effect of neighbourhood structure. These structures are represented in Fig 3 as (L5, L9, C9, C13, C21, and C25).
Transition rule. It is a set of rules that govern how each atom changes state. The next state of the current atom is decided by evaluating the status of neighbouring atoms. Let us take the best position of the neighbor M i around an atom as P t i ðM n Þ. The neighbour with the highest fitness value will be used to facilitate information diffusion and effectiveness. The transition rule can be defined as

PLOS ONE
where the present value of M i is represented by fitness À P t i ðM i Þ � . Therefore neighbour with the best fitness value from the CA structure can be obtained by using Eq (25). Later it will be used to update the state of the centre cell.

Particle swarm optimization
It is an evolutionary algorithm inspired by the behaviour of a bird's flock [31]. It offers a number of advantages, including a low mathematical complexity, a high capability for optimization, and easy implementation. PSO is comprised of two distinct processes. Birds carry out both processes. Initially, the bird will conduct a random search for the food source that is the closest to it. Later on, it will utilise its flying experience to determine the location of the meal. The shift in position has been referred to as velocity, and it varies with the passage of time. During the flight, particle speed increased stochastically towards its best point (personal best) and the community's best solution (global best) [32]. The candidate solution, which is represented by a particle, is a bird. In addition, food acts as a representation of the best possible solution to the problem. In this work i th particle is considered as .v id represent the velocity of each particle. With their initial velocity each particle start searching in the search space, where p i = p i1 , p i2 , p i3 . . ..p id represent the personal best position of each particle and p g = p g1 , p g2 , p g3 . . ..p gd represent the global best position of each particle. In PSO each particle updates its velocity and position by using the following equations Here c 1 and c 2 are acceleration coefficient which mostly controls the exploration and exploitation capability of the algorithm, inertia weight is represented by w, r 1 and r 2 , both represent the random numbers [0,1]. The fitness value is the determinant to analyse the quality of the best particle. The particle which has the best fitness value is taken as the global best solution [33].

Proposed hybrid algorithm
The framework of the proposed algorithm, as well as the improvement approach, have been discussed in detail in this section.
Framework of the MOCSHOPSO algorithm. The purpose of this research is to improve the performance of SHO by utilizing different techniques from CA and PSO for wellbore trajectory optimization problems. One crucial strategy is CA which is utilized to improve the performance of the SHO-PSO. The following three advantages have motivated to implement the hybridization strategy.
• The local search capability can be enhanced with the assistance of CA, as it enhances exploitation capability through interaction with its neighbours. However, the process of information transmission aids in exploration.
• As a result of the candidates' consequent solutions being attracted to the good SHO solutions, the SHO's convergence speed is quite fast, resulting in a local optima trap. The slow diffusion mechanism of CA in conjunction with SHO will aid SHO in avoiding the local optima trap.
• The velocity update mechanism of PSO will be utilized to improve the hunting mechanism of SHO. The velocity component of PSO is frequently managed by multiplying the particle's velocity by a factor. This regulation of velocity is intended to strike a balance between exploration and exploitation.
This programme generates semi-random populations of spotted hyenas. They are distributed in an n-dimensional lattice grid. The operational constraints are managed during the initialization of the agent's positions in the manner that has been proposed. The wellbore trajectory tuning variables x j are then initialised at random, but in a constrained environment, to get the desired results. If the initial positions satisfy the operational and non-negative constraints, they are accepted without further consideration. The values of the first tuning variable will be less chaotic and semi-randomly created as a result of this method. Then, in terms of wellbore trajectory optimization, the fitness value of each population is computed. Following that, it begins its update loop. It utilises the CA principle to create a new neighbourhood. During the neighbourhood generation process, some neighbours overlap. It enables the algorithm to incorporate an implicit migration mechanism. Additionally, it aids in the seamless diffusion of the best solutions throughout the population. As a result, it can retain a greater degree of diversity than the original SHO. Soft diffusion is critical in maintaining a balance between exploration and exploitation. The entire search space is divided into several sub-search regions by this approach. As a result, they can update the operation separately. However, if the neighbours overlap in this situation, information is transferred on an ad hoc basis.
However, the hunting mechanism of SHO has been modified by using the velocity update mechanism of PSO [34]. Therefore Eq (18) can be expressed as follows.
The updated hunting mechanism will be used by the proposed algorithm The pseudocode of the proposed algorithm has been expressed as follows. Adaptive neighbourhood. In the proposed method every search agent will create a neighbourhood according to the fitness value of the search agent. If the search agent has superior quality, then it will create a large size neighbourhood. It will increase the probability of overlapping with the best one. Besides, it will also increase the diversity of the solutions which also helps to avoid the local optima trap. This strategy is characterized as BL strategy [35]. On the other hand, the worst search agent will create a large neighbourhood to increase the chances of interacting with the best ones in the case of BS strategy [35]. It will increase the quality of the bad solutions. That means two approaches will provide the advantage to the algorithm either by avoiding the local optima trap or by improving the bad solutions. Besides both will maintain diversity. When the algorithm will try to create the neighbourhood for each search agent it will do so according to the following algorithm.

Wellbore trajectory optimization by MOCSHOPSO
Since wellbore trajectory optimization is a non-linear optimization problem, the spotted hyena populations have been initialized semi-randomly to make the approach less chaotic. After evaluating each search agent for the wellbore trajectory optimization problem, it selects the best search agents and clusters of all obtained solutions. Then it will create a neighbourhood by using an adaptive neighbourhood framework. It will also test the fixed structure like (L5, L9, C9, C13, C21, C25) for this problem. After creating a neighbourhood according to the CA concept, they are evaluated for optimization problems along with constraint violation checking. Later they exchange information among them. If any neighbour provides a better solution than the existing population (central cell) then it will update its position otherwise the central cell will be the best solution. This process will continue until the end criterion is met. To begin, a list of the six mentioned neighbourhoods is created. This list is sorted by neighbourhood size. The concept of radius is employed in this work to quantify the size of the neighbourhood [5]. The neighbourhoods in consideration are L5, C9, C9, C13, C21, and C25. Second, whenever the algorithm requests the neighbourhood of an individual, the adaptive approach selects the neighbourhood based on the quality of the provided individual.

Results & discussion
In this section, several statistical analyses are discussed to investigate the performance of the proposed algorithm qualitatively and quantitively for wellbore trajectory optimization. Additionally, the proposed hybrid method's performance for trajectory optimization was compared to that of the previously utilised metaheuristic algorithm MOCPSO [10]. To validate this work, real field data (secondary) from the Gulf of Suez oil field was used [4].

Comparative criterion
The evaluation of this work has been performed based on some statistical analysis. These parameters are the inverted generational distance (IGD), the spacing metric (SP), the maximum spread (MS), and the error ratio (ER). The next sections will describe these characteristics in detail.
Inverted Generational Distance (IGD). A researcher will obtain two types of Pareto solutions based on the solution obtained. The Pareto front solution is one type, whereas the approximate Pareto front solution is another. The term IGD is used to indicate this distinction between the two. However, a true Pareto front is necessary to calculate IGD. When the true Pareto front is unavailable, non-dominated solutions serve as a baseline against which to compare. The formulation of IGD can be derived as Here, the difference of Euclidean distance between the true Pareto front (Q) and approximate (Q � ) Pareto front is represented by d(Q, Q � ). It is a distance-based accuracy metric. If any algorithm gives a minimum value of IGD among the comparable algorithms, then it means it has obtained a Pareto front which is very near to the true Pareto front.
Spacing metric (SP). When the beginning and end of a Pareto front are unknown, it is required to determine the distribution of solutions. As such, SP is a metric that quantifies the distance variance between neighbouring vectors in the resulting non-dominated solutions. The minimum value of SP means the distance of solutions is comparatively less It indicates that solutions that are not dominated have a more balanced distribution. The formulation of the metric is as follows.
SP ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 1 jP À 1j Where, Maximum spread (MS). This metric is usually used to calculate the diversity of the obtained solutions and the coverage area of the solution. It is usually expressed as a difference between the two boundary solutions.

MS ¼
ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi Here, the maximum values in i th the objective is represented by a i and minimum values are represented by b i .
Error ratio (ER). This metric defines the number of solutions that are not secured a place in the obtained Pareto optimal set. It can be mathematically expressed as follows.
Where, the number of vectors in the currently available non dominated set is represented by n, e i = 0 if i� Pareto optimal set otherwise e i = 1. But ER = 0 is an ideal case that does not mean all solutions obtained by the proposed algorithm are within Pareto optimal set.

Comparative analysis
The statistical and analytical performances of the proposed algorithm on wellbore trajectory optimization are discussed in this section. In this analysis, the proposed algorithm is compared with the previously used algorithm MOCPSO and several other state-of-the-art algorithms such as MOSHO [21], MOCGWO [23]. Different algorithms utilized various parameters during optimization which are tabulated in Table 2. These algorithms tuned 17 variables within their constraint limit in order to produce the optimal trajectory.
Each of the algorithms under consideration obtained a Pareto front solution, as illustrated in Fig 4. As illustrated in the figure, MOCSHOPSO's Pareto front is more convergent and distributive than those of other algorithms. The blue circle depicts the evolved search agents, while the diamond shape (red colour) reflects the archive's non-dominated solutions. Because CA is hybridised, each search agent can communicate with its neighbours. It enhances the algorithm's potential for local search, which ultimately enhances the algorithm's convergence capability. Additionally, as illustrated in Fig 4, the proposed algorithm's Pareto front has a more diverse distribution of solutions than MOSHO (Fig 4C).
This substantial improvement is attributable to the addition of CA. By employing the characteristics of an adaptive neighbourhood, it raises the likelihood of better exploration, which ultimately aids the algorithm in providing a solution that is superior to the original. Additionally, CA's adaptive behaviour enhanced the likelihood of improving the worst solution. However, the quality of all derived Pareto fronts is statistically assessed in this work. Each method is run 30 times during this investigation, and all box plots are generated based on the collected data. The box plot more precisely expresses the distribution of solutions than any other style of plot. Additionally, it demonstrates the data sets' dispersion and symmetry [36]. Fig 5 illustrates the box plot for the IGD comparison metric. Whereas Table 3 contains the comparative data from the IGD analysis. The proposed approach achieved a mean value of 0.01283, which is the lowest value attained by any of the compared algorithms. Additionally, it outperformed the previously utilized MOCPSO algorithm for wellbore trajectory optimization. Furthermore, the data demonstrate that, althouth MOSHO has a larger mean value than MOCSHOPSO, but it's obtained best solution (0.00214) is better than MOCSHOPSO. The minimum IGD value also underscored the fact that the algorithm's Pareto front is closer to the true Pareto front than others. Due to the diffusion process of CA, it enabled the algorithm in exploring more search regions and gathering more isolated minima than was previously achievable. As a result, the gap between the obtained and true Pareto fronts has narrowed. Table 4 summarises the SP analysis's findings. The proposed technique achieved the lowest mean value among comparable algorithms (77.0527). This shows that the solutions generated by this algorithm are very evenly spaced. Throughout the hunting process, the recommended algorithm's search agents used PSO's velocity update feature to update their positions. This aided the algorithm in producing this positive outcome. The decrease in the value of h ! and F !

PLOS ONE
has also played a role in this problem. Fig 6 illustrates the box plot of the spacing metric analysis.
The proposed technique is also analysed in terms of MS and ER metrics in this work. The MS and ER experimental results are summarized in Tables 5 and 6, respectively. MS and ER box plots are presented in Figs 7 and 8, respectively.
As indicated by the experimental findings, MOCSHOPSO has the highest spread and the lowest error ratio among the comparison algorithms. The proposed algorithm achieves a mean value of 86.5454 for the MS, which is much greater than the value obtained by existing algorithms. This is due to the enhancement of MOSHO's exploration capability. Additionally, it was noted that the proposed algorithm more uniformly distributes non-dominated solutions than others. The proposed method results in a mean of 0.1122 for ER. This suggests in certain search regions, it has a lower proportion of significant non-dominated solutions than in

PLOS ONE
others. Additionally, it has been demonstrated that MOCSHOPSO is capable of collecting more isolated significant minima than other algorithms. This favourable finding also demonstrates the critical nature of including the slow diffusion mechanism of CA. Because an algorithm can achieve the lowest ER value only when the number of missing non-dominated solutions is minimized.
However, the proposed method has been compared qualitatively. Fig 9 depicts the data from each repository in terms of three objective functions. Before plotting, the values in the repository are sorted. MOCSHOPSO surpassed all other algorithms in terms of TMD, as shown in Fig 9(A). Though MOCPSO began with a lower value, MOCSHOPSO eventually reached the optimal minima faster.
From a similar analysis of torque, it is clear that the torque generated by MOCSHPSO is significantly less than that calculated by other methods. MOSHO performs similarly to the previously used algorithm MOCPSO in this case, while it performs better in repositories 10 and 20. Additionally, the suggested approach outperforms the comparative algorithm in terms of well profile energy. This suggests that the design parameter obtained will result in a less complex and cost-effective trajectory design.

PLOS ONE
Throughout the optimization process, the algorithms have taken into account the geological and operational constraints associated with wellbore trajectory optimization, as shown in Table 1. The algorithmic management of constraints and their required boundaries is depicted in Fig 10. It is shown that MOCSHPSO managed constraints effectively and with very little fluctuation. This is because the initialization is semi-random rather than random. It has resulted in less chaotic handling of constraints.
The performance of MOCSHOPSO with different neighbourhood structures. Six distinct variations of MOCSHOPSO are described in this section, each based on a different neighbourhood topology in addition to the adaptive one. These versions were evaluated against one another in terms of performance and statistical analysis. However, for wellbore trajectory optimization, MOCSHOPSO with C21 performed marginally better than its other versions in the case of fixed topology. As a result, C21 is recognized as a suitable neighbourhood structure for forming the optimal MOCSHOPSO. However, when adaptable neighbourhoods are used, they slightly outperform all fixed neighbourhood topologies. Table 7 summarizes the IGD metric analysis for each neighbourhood topology, and Fig 11 illustrates the box plot. According to  Table 7, the mean value for adaptive neighbourhoods is 0.3251 which is much less than the mean value for fixed neighbourhoods. C21 has a mean value of 0.5975 among the fixed topologies. As illustrated in Fig 11, adaptive neighbourhood (AN) has better distribution than others. Additionally, the Friedman test was used to determine the optimal neighbourhood topology among the neighbourhood topologies. It is a frequently used statistical test for multiple comparisons that is inherently non-parametric [37]. Table 8 has tabulated the data for the Friedman test analysis.

PLOS ONE
According to the Friedman test data, it can be concluded that AN performed better than other topologies. This is the optimal neighbourhood topology for optimizing the wellbore trajectory.
Radius also plays an important role during neighbourhood construction. If the radius becomes too small then it may fall into the local optima trap and may miss some significant minima. On the other hand, if it becomes too large it gives more attention to global search which may lead to slow convergence.
In this paper, various radius values are evaluated to determine the optimal radius for the neighbourhood topology. Table 9 contains the experimental data for IGD analysis, and Fig 12  depicts the box plot. According to the experimental results, the optimal radius for neighbourhood formation is R = 1.5. It has a mean of 0.2166, which is the smallest number among the compared radii.
Sensitivity analysis. It is a qualitative analytical technique used to ascertain the significance of tuning parameters utilized throughout the optimization process. Different tuning parameters contributed identically to attaining the minimal value for three optimization objectives throughout the optimization. There are seventeen tuning variables in this work, each with its constraints. This study will determine which parameter contributed the most and which one contributed the least. Among the two types of sensitivity analysis, local sensitivity analysis expresses the relevance of each variable but does not reveal any correlation information in the case of multi-objective optimization. However, with the global sensitivity analysis, adjustment of any variable, a correlation between all the objectives is recognized. Due to the nonlinear nature of wellbore trajectory optimization, global sensitivity analysis was performed in this case. In this paper, the Spearman correlation coefficient analysis was used to determine the global sensitivity [24]. Table 10 summarizes all experimental data. If a variable achieves the

Conclusion
In this work, a noble hybrid algorithm named, MOCSHOPSO has been proposed and demonstrated for wellbore trajectory optimization considering the TMD, well profile energy and torque. The hybridization is done with the incorporation of cellular automata (CA) and particle swarm optimization (PSO). The incorporation of CA has significantly improved the exploration capability of the algorithm. Besides the velocity update mechanism of PSO significantly contributed to the hunting mechanism of MOCSHOPSO. The performance of the proposed

PLOS ONE
algorithm in wellbore trajectory optimization is compared with the MOCPSO, MOSHO, and MOCGWO. This algorithm showed superior performances than all the three established algorithms. The significant improvements of the original algorithm have been validated through different statistical analyses. It has achieved the lowest IGD, SP and ER value on the other hand it has also given maximum value for MS. The hybridization has also significantly improved the capability of tracing the isolated minima which are validated by the lower value of ER. Moreover, an adaptive neighbourhood mechanism has proposed which has shown better performance than the fixed neighbourhood topologies. Additionally, the proposed algorithm has shown significant improvements during geophysical and operational constraints

PLOS ONE
handling of the wellbore trajectory optimization. Semi-random initialization of the search agents has helped the algorithm to achieve this less chaotic constraint handling. Moreover, this algorithm has provided optimum parameters for trajectory design which will help the industry to design a cost-effective and less complex trajectory design.  Author Contributions